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We have studied the cosmological evolution of domain wall networks in two, three and four spatial 
dimensions using high-resolution field theory simulations. The dynamical range and number of our 
simulations is larger than in previous works, but does not allow us to exclude previous hints of 
deviations to the naively expected scale-invariant evolution. These results therefore suggest that 
the approach of domain wall networks to linear scaling is a much slower process than that of cosmic 
strings, which has been previously characterized in detail. 



I. INTRODUCTION 

The cosmological consequences of primordial phase 
transitions that are thought to have happened in the 
early universe have been the subject of many studies in 
recent years. One of these inevitable consequences is the 
formation of topological defects P, |2| ■ While most such 
studies have focused on cosmic strings, domain walls can 
be of interest too, although observational constraints rule 
them out if their symmetry breaking scale is ry > IMeV 
0. However, note that this constraint is on the assump- 
tion that they are topologically stable and that their cos- 
mological evolution is scale- invariant {id est, that the 
network correlation length is proportional to the hori- 
zon size, or equivalently that the total area density is 
inversely proportional to it). 

In this work we want to highlight the fact that, in 
spite of some simple analytic arguments in favor of it 
U, this hypothetical scale-invariant solution has not yet 
received any strong backing from numerical simulations. 
In fact, in previous simulations of cosmological domain 
wall evolution 0, IE S IS 13 some hints for deviations 
to a scale-invariant evolution were found. The question 
we wish to address presently is whether these deviations 
remain if one performs larger runs. 

We report on a large set of high-resolution simulations 
of domain walls in two, three and four spatial dimen- 
sions, usinff the standard Press-Ryden-Spergel (PRS) al- 
gorithm p • Our boxes are typically larger (hence can be 
evolved for a longer dynamic range) that those in previ- 
ous works. Gains in statistical accuracy also come from 
performing large numbers of runs. We highlight the dif- 
ferences between running a large number of small runs, or 
a small number of larger ones. We do find some statistical 
evidence for deviations to the scale-invariant behaviour, 
though there is some dependency on dimensionality and 
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box sizes. While 2D and 3D numerical simulations have 
been previously carried out (though usually in smaller 
boxes or with smaller dynamic range) by a number of 
authors |E IE IE 1 we have also carried out for the 
first time a series of 4D simulations, which in a simple 
phenomenological way may be of interest to brane world 
scenarios 10', llj. 

In what follows we start with a brief description of the 
PRS algorithm and then proceed to present some tests of 
our code and the main results of our simulations. We also 
compare these with the work of previous authors, and 
highlight issues having to do with the slightly different 
methods used in each case. Finally, we briefly comment 
on the cosmological implications of our results. 

II. EVOLUTION OF THE DOMAIN WALLS 

Let us consider the evolution of a domain wall net- 
work in a flat homogeneous and isotropic Friedmann- 
Robertson- Walker (FRW) universe with line element 

ds^ = dt^ -a^{t){dx^ +dy^ + dz^), (1) 

where a(i), is the cosmological expansion factor, x, y 
and z are co-moving coordinates and t is physical time. 
The dynamics of a scalar field (p is determined by the 
Lagrangian density 

C = i0,a0^" - V{4,) , (2) 

where we will take F(0) to be the generic (fi^ potential 
with two degenerate minima given by 

ViP) = V.[^^^l)\ (3) 

This obviously admits domain wall solutions. By varying 
the action 

[ dt [ d^xy/^C , (4) 
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with respect to </> we obtain the field equation of motion: 

(5) 



52<^^„„a</. 2 dv 



where V is the Laplacian in physical coordinates and 
H = (da/dt)/a is the Hubble parameter. 

Following [a| we have modified the equations of motion 
in such a way that the co- moving thickness of the domain 
walls is fixed in co- moving coordinates. This has a small 
impact on the large scale dynamics of the domain walls 
and should not affect the quantities we want to measure 
for the purpose of testing scaling properties , provided 
a minimum acceptable thickness is used — this will be dis- 
cussed in Sect. IIIII On the other hand, it allows us to 
resolve the domain walls through the network evolution. 
For a discussion of analogous issues in the context of cos- 
mic strings see [T^ . 

Then equation |(SJ) becomes: 



din a 
dhirj J drj 



,dV 



(6) 



where ry is the conformal time and Pi and (32 are con- 
stants: /3 = is chosen in order to have constant co- 
moving thickness and a = 3 to require that the momen- 
tum conservation law of the wall evolution in an expand- 
ing universe is maintained 5] In fact we have numerically 
tested this with a — 2 and a — A and no significant dif- 
ferences were found in the quantities we want to measure 
in this context. 

Equation © was integrated using a standard finite- 
difference scheme, with the value of the Hubble damping 
term determined by numerical integration of the Fried- 
mann equation: 

d-q \ 1 

da Ho a^f7„oa-i + f^Aoa^ + (1 - ^^mo - ^Ao) 

(7) 

The conformal time evolution of the co-moving corre- 
lation length of the network = ^ was determined. 
Following 5] the co-moving area A was determined by 
finding the neighbouring grid points where </> has differ- 
ent signs and adding in these cases an increment AA di- 
vided by a weighting factor corresponding to the spatial 
gradient of 0: 



A 



OJ- 



V(f> 



d(j)/dx I + I d(t)/dy \ + \ dcjj/dz \ 



(8) 



where w = 1 if the link crosses the wall and a; = other- 
wise. Finally, the ratio between the kinetic and potential 
energy energy of the (j) field is roughly proportional to 
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E 



drj 



(9) 



This quantity is directly related to the root-mean squared 
velocity of the domain walls which is conserved if a scal- 
ing solution exists. We assume the initial value of cj) to be 



a random variable between —(f>o and +(f>o and the initial 
value of d(j)/dri to be zero. See Q for further discussion 
of these issues. 



III. SIMPLE TESTS OF THE CODE 

A total of several thousand simulations in two, three 
and four spatial dimensions were run for various box sizes 
and other initial conditions, as further detailed below. A 
number of observables were used as diagnostics for scal- 
ing. We start by describing the one that is closest to pre- 
vious work (alternatives will be described subsequently). 
In each run we looked for the best fit to the power law 
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where p^, is the wall energy density and rj the conformal 
time, and thereby measured the exponent A. In terms 
of the co-moving and physical correlation lengths of the 
network this corresponds to the following scaling laws, 
respectively 



1-A 



rj 



l-A/3 



(11) 



(12) 



Note that since we will be looking for deviations to the 
scale-invariant behaviour, which corresponds to A = 0, 
these will be larger if we work in co-moving coordinates, 
since any deviation exponent will be larger (though mea- 
surement errors could also differ). 

Each fit was made considering only the reliable dynam- 
ical range of the corresponding simulation. Note that in 
the beginning of the simulation the initial conditions will 
influence the evolution, and additionally the co-moving 
correlation length must be significantly larger than the 
wall thickness for the evolution to be sufficiently well de- 
fined. On the other hand, towards the end of the sim- 
ulation the boundary conditions may become relevant. 
Because in a (nearly) scale-invariant regime the average 
wall velocity is expected to be approximately constant, 
our choice of the reliable period was based on the con- 
stancy of the kinetic to potential energy ratio, as in . 

One finds that early in the simulations this ratio dis- 
plays strong oscillations while the initial conditions are 
relaxed. The relaxation timescale is typically given by 
the. wall thickness, which we parametrise by Wq The 
ratio then stabilizes and there follows a period where it is 
approximately constant. Finally, towards the end of the 
simulations it may start to decrease or increase slightly. 
Two examples are shown in Fig. ^ We do our fits in this 
intermediate, near-constant regime, selecting the appro- 
priate interval as follows. Firstly we aggressively choose 
a preliminary interval of conformal times, and determine 
the corresponding average value of the ratio. Secondly, if 
we find that at any point within the interval the kinetic 
energy varies by more than about ten percent, we reduce 
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FIG. 1: The evolution of the ratio between the kinetic and 
potential energy energy of the 4> field with conformal time 
for 2048^ simulations in the radiation (top) and matter (bot- 
tom) epochs. The simulations are run for a dynamical range 
of 1024. Wo is the wall thickness — notice that this sets the 
relaxation time scale. 



FIG. 2: The collapse of spherical domain walls with different 
thickness Wo- All simulations have a dynamical range of 1200. 
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FIG. 3: Measuring scaling exponents A for relatively small 2D 
(top set of points) and 3D (bottom set of points) simulations 
of domain wall networks. In each case 2^^ is the total number 
of points in a box, and the error bar is the standard deviation 
in an ensemble of 100 simulations. 



the interval (at either end, or both ends if necessary) un- 
til this is no longer so, and this final range is where the 
fit is then done. 

Let us pause at this point to comment on the influence 
of the wall thickness on the dynamics of the networks. 
Consider first the collapse of a single spherical domain 
wall of a given size, for different values of the thickness 
Wq. As we can see in Fig. |21 there is a minimum accept- 
able resolution of about 5 grid points (which will still 
yield inaccurate timescales), with 10 grid points being 
adequate. This can be confirmed by looking at network 
evolution, e.g. of 2D boxes in the radiation and matter 
epochs as in Fig. ^ Here the relevant figure of merit is 
the average ratio between the kinetic and potential ener- 



gies after the initial period of relaxation. 

We have carried out series of 100 test runs with differ- 
ent 2D and 3D boxes, to study the effect of the box size 
on the evolution of the networks. In each case we ran a 
number of different simulations, calculated an exponent 
for each simulation, and then computed the standard de- 
viation for A from the ensemble of boxes (the fit range 
used was the same for all simulations with the same dy- 
namical range). As shown in Fig. |21 the various results 
for the scaling exponent are consistent with each other, 
with the standard deviation in the result decreasing as 
we increase the box size. Also, as can be seen from the 
example in Fig. ^ the error function distribution for each 
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FIG. 4: The error function distribution for an ensemble of 100 
256^ matter era runs. The thin line is the Gaussian function 
distribution. 

ensemble of runs is approximately Gaussian (results for 
other series of runs are qualitatively similar) . Notice that 
the exponents differ in 2D and 3D — we shall return to 
this point in the main results section. 

Another simple test is to consider the evolution of the 
standard deviation in conformal time. In the scaling 
regime we would expect the standard deviation to be 
proportional to y/N/N cc rf'l'^^ N = {Lc/rff and Lc is 
the co-moving box size and 77 is the horizon size. If one 
allows for the the correction 7y~3/2(i-A)^ with A calcu- 
lated as before for all the simulation fit range, we found 
that the standard deviation was almost constant and in- 
dependent of conformal time. 



IV. NUMERICAL SIMULATION RESULTS 

Our results are summarized in Tabled where for com- 
parison purposes we also show a number of previously 
obtained results. In each case we ran a number of dif- 
ferent simulations, calculated an exponent for each sim- 
ulation, and then computed the standard deviation for 
A from the ensemble of boxes. Note that other authors 
calculate errors in slightly different ways. For example, 
first carries out an average of all runs and then does 
the fit on the 'averaged run'. 

We quote both the full dynamic range (in conformal 
time) of each simulation and also the dynamic range cor- 
responding to the interval that was actually used to fit for 
the scaling exponent. In some cases explicit details are 
not provided by previous authors and we have estimated 
them, either based on indirect information or by inspec- 
tion, e.g. of figures. For the case of the full dynamic 
range is effectively over-estimated, since early in their 
simulations dissipation was introduced in order to speed 



up the formation of the domain wall network. It is no- 
ticeable that larger, longer simulations do not necessarily 
correspond to a bigger dynamic range for the fit. Admit- 
tedly, this reflects the fact that different authors have dif- 
ferent criteria to decide when they trust their simulations 
to be near the scaling regime, and hence which range to 
use in the fits. By comparison with previous work, our 
method seems to be fairly conservative, and this is one 
of the reasons why we can get quite small error bars. We 
have verified through simple tests {e.g., changing the ten 
percent cutoff to three percent) that within error bars 
our results do not critically depend on it. 

Our results confirm previous work in that the approach 
to scaling is indeed very slow. We again find hints for 
deviations from the scale-invariant (A = 0) behaviour, 
which are consistent with most previous results (in par- 
ticular with the original runs of Press, Ryden and Spergel 
[^). This is to be contrasted this with the evolution of 
cosmic string networks, where field theory simulations 
with comparable or even smaller dynamic range find a 
relatively faster relaxation to scaling 0, On the 

other hand, notice that in the case of Garagounis and 
Hindmarsh no deviations from the scale-invariant be- 
haviour are found in their standard simulations, but a 
one-off very large 3D simulation does find them. We do 
agree with Q that there is not any strong evidence for 
a logarithmic correction to the scaling law (which has 
been suggested by 0), though it must be kept in mind 
that even with our larger dynamic ranges it would not 
be trivial to identify such a correction and, more to the 
point, to distinguish it from a power-law deviation to 
scale invariance. 

Moreover, we caution that comparisons between 
and the other simulations are not straightforward for var- 
ious reasons. Apart from their choice of physical rather 
than conformal time to present the results, they also use 
a different algorithm for calculating the domain wall ar- 
eas and, more importantly, dissipation is put into the 
early part of the simulations. While this has the obvious 
benefit of speeding up the formation of the wall network, 
and hence to some extent the approach to whatever is 
the late-time attractor solution for the network, it is not 
clear that it will not change the scaling properties them- 
selves. At a qualitative level, we would expect that his 
could be the difference between their results and the oth- 
ers. Also, for a generic discussion of the possible pitfalls 
of lattice-based simulations see . 

Since our simulations have a large dynamical range, 
one can make the point that having a single exponent to 
characterize the entire range might not be a very good 
approximation, for presumably if the network is indeed 
approaching a linear scaling regime then the deviations 
should be decreasing as the simulations evolve. This is 
indeed a valid objection, and a simple way of studying 
what happens is as follows. We have taken our sets of 
100 simulations of 4096^ and 256"^ boxes and divided the 
dynamic range of each into a number of equally spaced 
bins in conformal time (20 bins for the former, 12 for the 
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Dimension 


Reference 


Box Size 


Epoch 


»7o 


'7n 




Wo 


Full Range 


Fit Range 


Runs 


Scaling exponent A 


2D 


f5] 


1024^ 


Mat. 


1.0 


1.0 


1.0 


10 


512 


5.0 


10 


0.129 ±0.015 


2D 


f6] 


1024^ 


Rad. 


1.0 


1.0 


1.0 


5 


512 


10.0 


5 


0.12 ±0.04 


2D 


m 


1024^ 


Mat. 


1.0 


1.0 


1.0 


? 


512 


10.0 


few 


0.05 ±0.08 


2D 


[9] 


1024^ 


Rad. 


1.0 


0.1 


0.3 


3.3 


1536 


5.0 


5 


0.004 ±0.018 


2D 


f9] 


1024^ 


Mat. 


1.0 


0.1 


0.3 


3.3 


1536 


5.0 


5 


0.008 ±0.014 


2D 


This work 


4096^ 


Rad. 


1.0 


0.25 


1.0 


10 


2048 


200 


100 


0.0448 ± 0.0045 


2D 


This work 


4096^ 


Mat. 


1.0 


0.25 


1.0 


10 


2048 


200 


100 


0.0340 ± 0.0046 


2D 


This work 


8192^ 


Rad. 


1.0 


0.25 


1.0 


10 


4096 


500 


10 


0.0285 ± 0.0067 


2D 


This work 


8192^ 


Mat. 


1.0 


0.25 


1.0 


10 


4096 


500 


10 


0.0095 ± 0.0035 


3D 


f5] 


200^ 


Mat. 


1.0 


1.0 


1.0 


10 


100 


2.0 


9 


0.08 ±0.02 


3D 


f6] 


128^ 


Rad. 


1.0 


1.0 


1.0 


5 


64 


6.4 


5 


0.11 ±0.06 


3D 


m 


128^ 


Mat. 


1.0 


1.0 


1.0 


? 


64 


6.4 


few 


0.08 ±0.08 


3D 




512^ 


Rad. 


1.0 


0.1 


0.3 


3.3 


768 


5.0 


5 


0.006 ±0.013 


3D 


m 


512^ 


Mat. 


1.0 


0.1 


0.3 


3.3 


768 


5.0 


5 


0.003 ±0.012 


3D 


[9] 


4096^ 


? 


1.0 


0.08 


0.3 


3.3 


7680 


5.0 


1 


0.015 ± 0.003 


3D 


This work 


256^ 


Rad. 


1.0 


0.25 


1.0 


10 


128 


12.0 


100 


0.0811 ± 0.0044 


3D 


This work 


256^ 


Mat. 


1.0 


0.25 


1.0 


10 


128 


12.0 


100 


0.0681 ± 0.0032 


3D 


This work 


512^ 


Rad. 


1.0 


0.25 


1.0 


10 


256 


25.0 


10 


0.0343 ± 0.0023 


3D 


This work 


512^ 


Mat. 


1.0 


0.25 


1.0 


10 


256 


25.0 


10 


0.0468 ± 0.0019 


4D 


This work 


128* 


Rad. 


0.1 


0.05 


1.0 


10 


51 


5.0 


5 


0.0727 ±0.0013 


4D 


This work 


128* 


Mat. 


0.1 


0.05 


1.0 


10 


51 


5.0 


5 


0.0707 ± 0.0005 


4D 


This work 


128* 


Rad. 


1.0 


0.25 


1.0 


10 


64 


6.0 


8 


0.0743 ± 0.0025 


4D 


This work 


128* 


Mat. 


1.0 


0.25 


1.0 


10 


64 


6.0 


8 


0.0702 ± 0.0044 



TABLE I: The measured scaling exponent A for D = 2, D = 3 and D — 4 numerical simulations of domain wall networks. 
For comparison, we include both our own results and those of previous authors. The fourth column describes the cosmological 
epoch in which each set of simulations was performed. The ninth and tenth columns show respectively the total conformal 
time dynamic range of the simulation (or the range until the horizon reaches half the box size, for the case of simulations that 
extend beyond this) and the conformal time dynamic range of the part of the simulation that was actually used in order to fit 
for the scaling exponent. In some cases explicit numbers are not provided by authors; whenever possible we have estimated 
them. 



latter). We have then calculated a scaling exponent for 
each bin and run, and finally averaged over the ensemble 
of 100 runs for each of the bins. 

The results are shown in Figs. |31 and O A notice- 
able feature is that the error bars tend to increase as the 
dynamic range increases. This is because as the Hubble 
volume becomes a larger fraction of the box size there are 
effectively less independent regions to sample from. The 
plots of the 3D simulations do show some evidence of evo- 
lution of the scaling exponent, which (if one neglects the 
initial bin) starts out around Xi = 0.1 and then slowly but 
steadily decreases. The situation is not so simple in the 
2D case. Here in fact the exponent starts out at approx- 
imately the same value — a consequence of the fact that 
the same algorithm is being used to generate the initial 
conditions. However, in the longer dynamic range avail- 
able the exponent oscillates and eventually 'overshoots' 
the linear scaling regime, reaching negative exponents to- 
wards the end of the runs. Whether this is a numerical 
artifact or has some physical significance is presently un- 
clear. 



Notice that the fact that the exponent is negative (that 
is, ^ is growing proportionally faster than time) does not 
necessarily mean there is any violation of causality. By 
this time there can be a sufhciently large number of walls 
inside a Hubble volume for such a regime to ensue, albeit 
transiently. A similar situation happens with the Kib- 
ble regime for the evolution of cosmic string networks 
[TsllTa l. There, the network must be friction-dominated 
for such a regime to occur. In the present context this is 
not the case, though it can't be excluded that with the in- 
creasing dynamic range a background of small-scale fea- 
tures accumulates which could slow the walls down. Note 
that measuring velocities in field theory simulations is no- 
toriously difficult — see ^3 fo^' ^ discussion of the issue 
in the context of cosmic string simulations. 

Finally, it must be emphasized that if the network does 
indeed not have a scale-invariant solution, or if it does 
have one but it takes much longer to reach it than the 
dynamic range available to existing simulations, then the 
initial conditions with which one start the simulations 
will of course play an important role in determining the 
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FIG. 5: The evolution of the scahng exponent A with con- 
formal time for 100 4096^ simulations in the radiation (top) 
and matter (bottom) epochs. The dynamical range of each 
simulation was divided into 20 equally spaced bins, and an ex- 
ponent measured for each bin and then averaged over the 100 
runs. The error bar is the standard deviation in the ensemble. 

subsequent evolution of the networks, and may therefore 
affect the scaling exponents obtained. This also needs 
to be taken into account when comparing the results 
obtained by different authors. Understanding the spe- 
cific differences between the different available numerical 
pipelines, and possibly even carrying out independent 
tests with other algorithms are clearly key issues, and 
deserve further investigation. 



V. DISCUSSION AND CONCLUSIONS 

We have presented the results of high-resolution, long 
dynamic range field theory simulations of domain wall 

networks. In addition to the usual 2D and 3D numerical 
simulations, we have also performed some simulations of 
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FIG. 6: The evolution of the scaling exponent A with con- 
formal time for 100 256^ simulations in the radiation (top) 
and matter (bottom) epochs. The dynamical range of each 
simulation was divided into 12 equally spaced bins, and an ex- 
ponent measured for each bin and then averaged over the 100 
runs. The error bar is the standard deviation in the ensemble. 



domain wall evolution in four spatial dimensions. These 
can be relevant, in a phenomenological way, for models 
with additional space-time dimensions, such as the so- 
called brane world scenarios. 

Our results are further evidence of the fact that these 
networks have a very slow relaxation towards the ex- 
pected asymptotic linear scaling solution. Even with the 
longest dynamic ranges we can presently run, we can not 
find unambiguous evidence for the onset of a linear scal- 
ing solution. The large number and size of our simula- 
tions has allowed us to measure scaling exponents with 
considerably smaller error bars than in previous works, 
thereby strengthening the evidence for small violations 
of scale invariance in the ranges we can probe. However, 
one is always limited by the finite dynamical range of our 



7 



simulations and it is always possible that future simula- 
tions with a much larger dynamical range will produce 
results consistent with scaling. It is also worth emphasiz- 
ing that the measured deviations tend to decrease when 
we increase the box size. Box size effects, as well as 
possible dependencies on the way initial conditions are 
generated, clearly deserve further investigation. 

It is interesting to try to understand what could be the 
physical mechanism behind the violations to the scale 
invariance. Notice that typically they are such that 
the field configurations are equilibrating somewhat more 
slowly than would be allowed by causality (with the ex- 
ception discussed above). As was already pointed out, 
measuring velocities is a very difficult task in this type 
of simulation, but we do find that our networks are rel- 
ativistic. A possible explanation for the slower evolution 
is that it is due to the existence of long-range forces be- 
tween the walls. However, if that was the case we would 
expect that the exponents A, measuring the deviations 
from scale invariance, would decrease as the number of 
spatial dimensions increases: the more dimensions one 
has the easier it is for the walls to unwind, and the smaller 
(in relative terms) will be the effect of these forces. In 
fact, the opposite is the case in our simulations, although 
this is compounded by the fact that the dynamic ranges 
we can run are obviously smaller for larger dimensions, 
so one is not really comparing like with like. 

Finally, let us point out that the existence of viola- 
tions to scale invariance can have some cosmological con- 
sequences. Consider a network formed as early as ob- 



servationally allowed, that is (roughly speaking, for the 
simplest standard scenarios) when the temperature of the 
Universe (photon temperature) was T ^ O.lMeV. At for- 
mation there would be on average one domain per horizon 
volume, and that would also be the case today if the net- 
work quickly reached the linear scaling regime. However, 
if we assume that from then on = 'rj^~^ with A then 
the number of walls per horizon volume today could be 
as high as about ten for allowed values of the exponent. 
In other worlds, each Hubble volume could be divided 
into a significant number of domains, each conceivably 
having different properties. An example of such a model 
is described in jl7i|. Although these would be important 
modifications to the standard domain wall scenario we 
do not expect the limits on the wall mass per unit area 
(or equivalently to the symmetry breaking scale) to be 
severely modified. 
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